// Some code to esimate find the correlation function for
// equilateral triangles.  I used it for debugging purposes,
// but here it is in case anyone wants to use (or modify) it.

vector<double> sumepepep(nrbins,0.);
vector<double> sumepepex(nrbins,0.);
vector<double> sumepexex(nrbins,0.);
vector<double> sumexexex(nrbins,0.);
vector<double> sumwt(nrbins,0.);

void EqTri_4(vector<ThreePtData>& eqtridata, const Cell& c1, const Cell& c2, const Cell& c3, double d3)
{
  if (d3 < minsep) return;
  if (d3 > maxsep) return;

  double d2 = fabs(c1.MeanPos()-c3.MeanPos());
  if (d2-c3.Size() > 1.25*d3) return;
  if (d2+c3.Size() < 0.8*d3) return;

  double d1 = fabs(c3.MeanPos()-c2.MeanPos());
  if (d1-c3.Size() > 1.25*d3) return;
  if (d1+c3.Size() < 0.8*d3) return;

  if (c3.Size() > 0.1*d3) {
    EqTri_4(eqtridata,c1,c2,c3.Left(),d3);
    EqTri_4(eqtridata,c1,c2,c3.Right(),d3);
  } else {
    DComplex centroid = (c1.MeanPos() + c2.MeanPos() + c3.MeanPos())/3.;
    DComplex q1 = centroid - c1.MeanPos();
    DComplex q2 = centroid - c2.MeanPos();
    DComplex q3 = centroid - c3.MeanPos();
    double www = c1.Weight()*c2.Weight()*c3.Weight();
    const DComplex e1 = c1.WE()*SQR(conj(q1))/norm(q1);
    const DComplex e2 = c2.WE()*SQR(conj(q2))/norm(q2);
    const DComplex e3 = c3.WE()*SQR(conj(q3))/norm(q3);
    const DComplex e2e3 = e2*e3;
    const DComplex e2e3c = e2*conj(e3);
    const DComplex eee = e1*e2e3;
    const DComplex ecee = conj(e1)*e2e3;
    const DComplex eece = e1*conj(e2e3c);
    const DComplex eeec = e1*e2e3c;

    int k = int(floor((log(d3)-logminsep)/binsize));
    Assert(k>=0);
    Assert(k<nrbins);

    sumepepep[k] += e1.real()*e2.real()*e3.real();
    sumepepex[k] += e1.imag()*e2.real()*e3.real()/3.;
    sumepepex[k] += e1.real()*e2.imag()*e3.real()/3.;
    sumepepex[k] += e1.real()*e2.real()*e3.imag()/3.;
    sumepexex[k] += e1.imag()*e2.imag()*e3.real()/3.;
    sumepexex[k] += e1.imag()*e2.real()*e3.imag()/3.;
    sumepexex[k] += e1.real()*e2.imag()*e3.imag()/3.;
    sumexexex[k] += e1.imag()*e2.imag()*e3.imag();
    sumwt[k] += www;
    eqtridata[k].gam0 += eee;
    eqtridata[k].gam1 += ecee;
    eqtridata[k].gam2 += eece;
    eqtridata[k].gam3 += eeec;
    eqtridata[k].weight += www;
  }
}

void EqTri_3(vector<ThreePtData>& eqtridata, const Cell& c1, const Cell& c2, const Cell& c3)
{
  double d3 = fabs(c2.MeanPos()-c1.MeanPos());
  if (d3 + c2.Size() + c1.Size() < minsep) return;
  if (d3 - c2.Size() - c1.Size() > maxsep) return;

  if (c2.Size() > 0.1*d3) {
    if (c1.Size() > 0.1*d3) {
      EqTri_3(eqtridata,c1.Left(),c2.Left(),c3);
      EqTri_3(eqtridata,c1.Right(),c2.Left(),c3);
      EqTri_3(eqtridata,c1.Left(),c2.Right(),c3);
      EqTri_3(eqtridata,c1.Right(),c2.Right(),c3);
    } else {
      EqTri_3(eqtridata,c1,c2.Left(),c3);
      EqTri_3(eqtridata,c1,c2.Right(),c3);
    }
  } else {
    if (c1.Size() > 0.1*d3) {
      EqTri_3(eqtridata,c1.Left(),c2,c3);
      EqTri_3(eqtridata,c1.Right(),c2,c3);
    } else {
      EqTri_4(eqtridata,c1,c2,c3,d3);
    }
  }
}

void EqTri_2(vector<ThreePtData>& eqtridata, const Cell& c12, const Cell& c3) 
{
  if (c12.Size() < halfminsep) return;

  EqTri_2(eqtridata,c12.Left(),c3);
  EqTri_2(eqtridata,c12.Right(),c3);
  EqTri_3(eqtridata,c12.Left(),c12.Right(),c3);
}

void EqTri_1(vector<ThreePtData>& eqtridata, const Cell& c123) 
{
  if (c123.Size() < halfminsep) return;

  if (c123.Size() > outputsize) {
    recursen++;
    dbg<< string(recursen,',')<<"EqTri_1: size = "<<c123.Size()<<endl;
  }

  EqTri_1(eqtridata,c123.Left());
  EqTri_1(eqtridata,c123.Right());
  EqTri_2(eqtridata,c123.Left(),c123.Right());
  EqTri_2(eqtridata,c123.Right(),c123.Left());

  if (c123.Size() > outputsize) {
    dbg<<string(recursen,',')<<"Done EqTri_1\n";
    recursen--;
  }
}

// This section goes into main():

  vector<ThreePtData> eqtridata(nrbins);
  EqTri_1(eqtridata,*wholefield);

  ofstream eqtriout("eqtri.out");
  eqtriout <<"gam0 gam1 gam2 gam3\n";
  for(int i=0;i<nrbins;i++) {
    eqtriout << fix(minsep*exp((i+0.5)*binsize)/60.);
    eqtriout << sci(real(eqtridata[i].gam0)/eqtridata[i].weight);
    eqtriout << sci(imag(eqtridata[i].gam0)/eqtridata[i].weight);
    eqtriout << sci(real(eqtridata[i].gam1)/eqtridata[i].weight);
    eqtriout << sci(imag(eqtridata[i].gam1)/eqtridata[i].weight);
    eqtriout << sci(real(eqtridata[i].gam2)/eqtridata[i].weight);
    eqtriout << sci(imag(eqtridata[i].gam2)/eqtridata[i].weight);
    eqtriout << sci(real(eqtridata[i].gam3)/eqtridata[i].weight);
    eqtriout << sci(imag(eqtridata[i].gam3)/eqtridata[i].weight);
    eqtriout << sci(eqtridata[i].weight)<<endl;
  }
  eqtriout <<"+++  ++x   +xx   xxx\n";
  for(int i=0;i<nrbins;i++) {
    eqtriout << fix(minsep*exp((i+0.5)*binsize)/60.);
    eqtriout << sci(sumepepep[i]/sumwt[i])<<sci(sumepepex[i]/sumwt[i]);
    eqtriout << sci(sumepexex[i]/sumwt[i])<<sci(sumexexex[i]/sumwt[i]);
    eqtriout << sci(sumwt[i])<<endl;
  }

  int j=nubins-1;
  int k2=nvbins/2;
  int k1=k2-1;
  for(int i=0;i<nrbins;i++) for(int k=k1;k<=k2;k++) {
    const ThreePtData& bindata = data[i][j][k];
    if (bindata.weight > 0.) {
      eqtriout << fix(bindata.meanr/60.);
      eqtriout << fix(bindata.meanu);
      eqtriout << fix(bindata.meanv);
      double sig = sqrt(bindata.vargam);
      eqtriout << sci(real(bindata.gam0));
      eqtriout << sci(imag(bindata.gam0));
      eqtriout << sci(real(bindata.gam1));
      eqtriout << sci(imag(bindata.gam1));
      eqtriout << sci(real(bindata.gam2));
      eqtriout << sci(imag(bindata.gam2));
      eqtriout << sci(real(bindata.gam3));
      eqtriout << sci(imag(bindata.gam3));
      eqtriout << sci(bindata.weight) << endl;
    }
  }

